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Abstract 

We present a new calculation, in the nonrelativistic QCD (NRQCD) factorization formalism, 
of tlie relativistic corrections to the double-charmonium cross section a[e^e~ J/ip + r/c] at the 
energy of the Belle and BABAR experiments. In comparison with previous work, our calculation 
contains several refinements. These include the use of the improved results for the nonperturbative 
NRQCD matrix elements, the resummation of a class of relativistic corrections, the use of the 
vector-meson-dominance method to calculate the fragmentation contribution to the pure QED 
amplitude, the inclusion of the effects of the running of a, and the inclusion of the contribution 
that arises from the interference between the relativistic corrections and the corrections of next- 
to-leading order in a^. We also present a detailed estimate of the theoretical uncertainty. We 
conclude that the discrepancy between the theoretical prediction for a[e'^e~' J ftp + rjc] and the 
experimental measurements has been resolved. 
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I. INTRODUCTION 



For a number of years, one of the largest discrepancies in the standard model has been 
the disagreement between theory and experiment for exclusive double-charmonium process 
e+e~ — > J/t/j + rjc at the S-factory energy of 10.58 GeV. Initially, the Belle Collaboration 
reported for the cross section times the branching fraction into four or more charged tracks 
a[e^e~ — > J /if) + r]!^ x 5>4 = 33lg ± 9 fb (Ref. [1]). The first theoretical predictions were 
based on NRQCD factorization calculations [2] at leading order in ccs, the QCD coupling 
constant, and the heavy-quark (or antiquark) velocity in the quarkonium rest frame. 
These predictions were (T[e"^e~ — > J/'4> + 77c] — 3.78 ± 1.26 fb (Ref. [3]) and a[e^e~ — > 
J/^ + 77c] = 5.5 fb (Ref. [4]).^ The calculation of Ref. [3] includes QED effects, while that of 
Ref. [4] does not. Other differences between these calculations arise from different choices 
of the charm-quark mass mc, NRQCD matrix elements, and a^. The sensitivities of the 
calculations to these choices are indicative of large sources of uncertainty in the theoretical 
calculations that have not yet been quantified. 

More recently, the Belle Collaboration has measured the production cross section times 
the branching fraction into more than two charged tracks and finds that (7[e^e~ — > J/'^ + 
ryj X B>2 = 25.6 ±2.8 ±3.4 fb (Ref. [5]). The BABAR Collaboration has also measured this 
quantity, and obtains a[e~^e~ J/ip + rjc] x = 17.6 ±2.8 ±2.1 fb (Ref. [6]). These new 
experimental results have narrowed the gap between theory and experiment. 

An important recent theoretical development is the calculation of the corrections of next- 
to- leading order (NLO) in (Ref. [7]). These yield a K factor of about 1.96. While this 
K factor is substantial, it does not, by itself, completely remove the discrepancy between 
theory and experiment. 

Relativistic corrections to a"[e+e~ — > J/ip + rjc] also make a significant contribution to the 
theoretical prediction. These corrections arise in two ways. First, they appear directly in 
the corrections of order and higher to the process e+e~ — > J/ip + r]c itself. Second, they 
arise indirectly when one makes use of phenomenological determinations of certain NRQCD 
matrix elements that appear in the expression for a[e'^e~ — > J/ip+rjc]- For example, the rele- 
vant matrix element of leading order in v for the J/ ip can be determined phenomenologically 

^ The authors of Ref. [3] initially reported a cross section of 2.31 ± 1.09 fb, but later corrected a sign error 
in the QED interference term to arrive at the value cited above. 
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from the experimental value for the width for J/%1) ^ e'^e~ and the theoretical expression 
for that process. However, the theoretical expression contains relativistic corrections, which 
then indirectly affect the calculation of a[e^e~ J/ip + r]^. The first relativistic correction 
appears at order v"^. (f^ 0.3 for charmonium.) In Ref. [3], the order-f^ correction was 
calculated and was found to be about 1.95(w^)j/^ + 2.37(f^)^^. Here, {v'^)h is the ratio of 
an order-v^ nonperturbative NRQCD matrix element to the leading-order matrix element 
in the quarkonium state H. The large coefficients in the order-T;^ correction potentially lead 
to a relativistic correction. In Ref. [3] , the K factor for the relativistic corrections was found 
to be 2.0li'\^. The large uncertainties arose from large uncertainties in the NRQCD matrix 
elements. 

Recently, progress has been made in reducing the uncertainties in the order-u^ NRQCD 
matrix elements by making use of a potential model to calculate the quarkonium wave 
function [8]. The results of Ref. [8] allow one to make a meaningful prediction for the 
relativistic corrections to a[e^e~ — > J/-^ + 77c]- Making use of these results to compute the 
relativistic corrections and taking into account the corrections of NLO in cts, the authors of 
Ref. [9] have given the prediction a[e^e~ — >• J/V' + r]^ = 17.5 ± 5.7 fb. 

The authors of Ref. [10] have taken a different approach, determining the NRQCD matrix 
elements of leading order in v and of relative order v'^ by using T[J/%1) —>■ e+e"], r[r]c 77], 
and r[J/'0 — > hght hadrons] as inputs. Their result, a[e^e~ — > J/ip + rjc] — 20.04 fb, is 
in agreement with the result of Ref. [9]. However, as we shall discuss, the values of the 
individual matrix elements that were used in Ref. [10] differ significantly from the values 
that were used in Ref. [9] . 

The results of Refs. [9] and [10] suggest that there is no longer a discrepancy between 
the experimental measurements and the theoretical prediction. Nevertheless, it is useful to 
include further refinements that improve the precision of the theoretical prediction and to 
estimate as precisely as possible the various theoretical uncertainties. 

In the present paper, we carry out a new calculation of the relativistic corrections to 
(T[e+e~ J/tp+rjc]. We include the effects of pure QED processes, as well as QCD processes. 
In the case of the pure QED processes, we incorporate a further refinement by making use 
of the vector-meson-dominance (VMD) formalism to compute the photon-fragmentation 
contribution. This approach reduces the theoretical uncertainties that are associated with 
the pure QED contribution. In our calculation, we make use of the approach of Ref. [8] to 
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resum a class of relativistic corrections to all orders in v. We also compute the contribution 
that arises from the interference between the relativistic corrections and the corrections of 
NLO in as- Our calculation takes advantage of the new higher-precision determinations of 
the relevant NRQCD matrix elements in Ref. [11]. We make use of the detailed error analysis 
of Ref. [11] to estimate the theoretical uncertainties in our calculation, some of which are 
highly correlated. 

The remainder of this paper is organized as follows. In Sec. II, we discuss the general 
form of the amplitude for the process e'^e~ — > J/ip+rjc, and, in Sec. Ill, we discuss the corre- 
sponding quark-level amplitude. Sec. IV contains the expression for the NRQCD expansion 
of the amplitude and a discussion of the matching between NRQCD and full QCD. In Sec. V, 
we describe the rcsummation method that we employ. Sec. VI contains the specifics of the 
frame and coordinate system that we use in the calculation. We present explicit formulas 
for the cross section in Sec. VII. The VMD method for computing the fragmentation contri- 
bution to the pure QED amplitude is summarized in Sec. VIII. We specify how we choose 
quarkonium masses in the calculation in Sec. IX. We present the method that we use to 
compute the interference between the relativistic corrections and the corrections of NLO in 
as in Sec. X. We give our numerical results in Sec. XI and compare them with the results 
from previous calculations in Sec. XII. Finally, we summarize and discuss our results in 
Sec. XIII. 



II. GENERAL FORM OF THE AMPLITUDE FOR e+e" J/^ + rj^ 

Let us consider the amplitude for the exclusive process 7* J/ip{Pi, A) + ric{P2), where 
A is the helicity of the J/ip. The ^'-matrix element for e+(/c2)e"(/ci) —>■ J/ip{Pi, A) -f- rjc{P2) 
can be written as 

M{\)^L''A^[J/'4^{X)+r^ci (1) 
where the leptonic factor is defined by 

L^^ ^ -i^v(k2)ru{h). (2) 
s 

Here, Cc is the electric charge of the charm quark and s — 4:El^^^ is the square of the e'^e" 
center-of-momentum (CM) energy. ^^[J/'0(A) -|- rjc] in Eq. (1) is the vacuum-to- J/'^ -|- rjc 
matrix element. It can be expressed in the following form, which derives from the Lorentz 



FIG. 1: Feynman diagrams for the process e~^e~ — > J/tp + rjf. at leading order in a^. The wavy 
line represents a photon, the curly line represents a gluon, and the straight lines represent the 
leptons and heavy quarks. There are six additional diagrams that can be obtained by reversing 
the directions of the arrows on the heavy-quark lines and/or by replacing the gluon by a photon. 

invariance of the amplitude and the parity conservation of the strong and electromagnetic 
interactions [3]: 

A^[J/^{X) + Vc] = {J/^{Pi, A) + Vc{P2)\U0m = tAe,,^pP^P^e*''{X), (3) 

where J^(0) is the electromagnetic current, e*(A) is the polarization four-vector of the J/ip 
with helicity A whose components in the J/ip rest frame are e*(A) = [0,e*(A)]. The conven- 
tion for the antisymmetric tensor in Eq. (3) is chosen so that 60123 = +1- We note that A is 
parity invariant and that A'^ transforms as a four-vector under parity. 

III. AMPLITUDE FOR 7* QQi^Si) + QQCSq) 

The exclusive process 7* — > J/ip{Pi,X) + rjc{P2) involves the decay of a virtual photon 
into two heavy quark- ant iquark {QQ) pairs Q{Pi)Q{pi) {i = I or 2). Both of the pairs are 
in color-singlet states. At leading order in as, the process proceeds through the diagrams 
shown in Fig. 1, plus two additional diagrams in which the directions of the arrows on the 
heavy-quark lines are reversed. There are also purely electromagnetic contributions to the 
process. At leading order in the QED coupling a, two types of diagrams contribute to the 
QED processes. The first type consists of the diagrams shown in Fig. 1 (plus two others in 
which the directions of the arrows on the heavy-quark lines are reversed) , but with the gluon 
replaced by a photon. The second type consists of diagrams in which a photon fragments 
into a J/ip- One of these diagrams is shown in Fig. 2. (There is an additional diagram in 
which the directions of the arrows on the heavy-quark line on the r]c side are reversed.) 
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FIG. 2: Feynman diagram for the process e^e^ — > J /ip + r/^ in which a photon fragments into a 
JZ-i/^. There is an additional diagram that can be obtained by reversing the directions of the arrows 
on the heavy-quark hne on the r^c side. 

A. Kinematics 

The Q and Q that evolve into the charmonium Hi with momentum Pi have momenta pi 
and pi, where 2 = 1 denotes the J/ip and i = 2 denotes the rjc- The Q{pi)Q{pi) pair is in 
a spin-triplet S'-wave state and the Q{p2)Qij>2) pair is in a spin-singlet S'-wave state. The 
four-momenta of the Q and Q in the z-th pair are expressed in terms of the total momentum 
Pi and the relative momentum gj: 

Pi = \Pi + gi, (4a) 
Pi = \Pi - qi. (4b) 

Pi and Qi are chosen to be orthogonal: Pi ■ qi = 0. In the rest frame of the i-th QQ pair, 
the explicit components of the momenta listed above are Pi = [2E{qi),0], qi = {0,qi), 
Pi = [E{qi),qi], and pi = [E{qi), -Qi], respectively, where E{qi) = ^/ml + qf is the energy 
of the Q or the Q in the QQ rest frame. 

B. Spin and color projectors 

A production amplitude of a Q{pi)Q{pi) pair can be expressed in the form 

u{pi)Av{pi) = Tr [A v{p,i)u{pi)] , (5) 

where ^ is a matrix that acts on spinors with both Dirac and color indices. The amplitude 
in Eq. (5) can be projected into a particular spin and color channel by replacing v{pi)u{pi) 



6 



with a projection matrix. The color projector tti onto a color-singlet state is 

= (6) 



where 1 is the 3x3 unit matrix of the fundamental representation of SU(3). The color- 
singlet projector (6) is normalized so that Tr[7ri7r|]=l. The projector of the pair Q{pi)Q{pi) 
onto a spin-triplet state with helicity A and the projector of the pair Qip2)Q{p2) onto a 
spin-singlet state are denoted by n3(pi,pi, A) and Ili{p2,p2), respectively. The projectors, 
valid to all orders in q^, are given in Ref. [12]: 

n3(pi,Pi,A) = - I -(p,-m,)^*{Xm + 2E{q,)]{^, + m,), (7a) 

MP2,P2) = ^ . , -Ap2 - mch'[r2+2E{q2W2 + me), (7b) 

AV2E{q2)[E{q2) +mc\ 

where the spin-polarization vector e*(A) satisfies Pi-e*(A) = 0. The spin projectors in Eq. (7) 
are normalized so that 

Tr[n3(pi,pi, A)n^(pi,pi, A)] = 4p?p?, (8a) 

iv[ni(p2,P2)nl(p2,P2)] = 4p^p^. (8b) 

Since we are considering an exclusive process, in which no hadrons are present other than 
the J/ip and the rjc, we consider only the states of the QQ pairs that have the same quantum 
numbers as the J/ip and the rjc- That is, the pair Q{pi)Q{pi) must be in a color-singlet, 

spin-triplet iS-wave state, as is the J/ip, and the pair Q{p2)Q{p2) must be in a color- singlet, 
spin-singlet ^-wave state, as is the ric- 

The spin projectors in Eq. (7) can be simplified as follows: 

ni(p2,P2) = 2^E{q2) ~ "^"^ ^ ' ^^^^ 

The spin projectors in Eq. (9) are also valid to all orders in qi. 

In addition, we provide the following formulas, valid to all order in q^, which are useful 
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in this calculation: 



7an3(pi,Pi, A)7° 



rricypi - pi ) 
2[E{qi)+mc 



rrir 



(10a) 
(10b) 



C. Projections of the four-quark states 

Prom the full QCD amphtude for 7* Q{pi)Q{Pi)Q{p2)Q{p2) , one can project out the 
amplitude for 7* — > QQ{Pi,qi,X) + QQ{P2,q2), where the first pair is in a color-singlet, 
spin-triplet state with helicity A and the second pair is in a color-singlet, spin-singlet state. 
Applying the spin projections to both QQ pairs simultaneously, one obtains 

^^(Pi, gi. A; P2, q2) = TV{^'^[7* ^ Q{pi)Q{pi)Q{p2)Q{p2)\ [Mpi.Pi. \)mi\®[MP2.p2)mi\] , 

(11) 

where ^^^[7* ^ Q{pi)Q{pi)Q{P2)Q{p2)] is the full QCD amplitude for 7* ^ 
Q{pi)Q{pi)Q{p2)Q{p2) , Ai is the vector index of the virtual photon, and Aq{Pii qi, A; P2, ^2) 
is the amplitude for 7* — >• QQ{Pi, qi, A) + QQ{P2, 52)- 

In the amplitude (11), the QQ pairs are not necessarily in S'-wave orbital-angular- 
momentum states. One can project out the S'-wave amplitude by averaging, for each QQ 
pair, over the direction of the relative momentum qi in the QQ{Pi) rest frame. The amplitude 
for 7* ^ QQi^S^, Pi, A) + QQi^So, P2) is 

^^(^^1, Pi, A; ^^0, P2) = ^^(Pi,gi,A;P2,g2), (12) 

where the bar on the right side of Eq. (12) is the average over the angles of both qi and q2 
in the Pi and P2 rest frames, respectively: 

^^(Pi,gi,A;P2,g2)= j ^|^^^(Pi, gi. A; P2, 52). (13) 

(Kti is the solid-angle element of q^, defined in the Pj rest frame. Once we have averaged 
over angles, the dependence in the amplitude (12) reduces a dependence only on q\ and 
qf|. Note that Pj depends on implicitly: Pf — A{rn?^ + qf)- 
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IV. NRQCD EXPANSION OF THE AMPLITUDE AND MATCHING 



The NRQCD expansion of Eq. (3) in terms of the vacuum-to- J/ and vacuum-to-77c 
matrix elements is 

^'^[J/V'IA) +77c] = V2^V2^5]ci^;,„(J/V'(A)|C»^|0)(77c|C»n|0), (14) 

m,n 

where the are short-distance coefficients and the Om and the 0„ are NRQCD operators. 
The quantities mi and m2 represent the J/%1) and r]c masses, respectively. However, as we 
shall discuss later, these quantities are not necessarily equal to the physical meson masses, 
but may instead be expressed as functions of the heavy-quark masses via the nonrelativistic 
expansion of NRQCD. The factor \j2mi\/2m2 appears on the right side of Eq. (14) because 
we use relativistic normalization for the meson states in A'^[J / 'il){X) -\- r^J, but we use con- 
ventional nonrelativistic normalization for the NRQCD matrix elements on the right side of 
Eq. (14). 

Now we approximate the formula (14) by retaining only those operator matrix elements 
that connect the vacuum to the color-singlet QQ Fock states of the quarkonia. Then, we 
have 

oo oo 

^^[J/^(A) + = v^V2^5]^c^JA)(J/^(A)|^t(_iB)2™a-e(A)x|0) 

m=0 n=0 

X(,7,|V'^(-|B)2"X|0), (15) 

where the short-distance coefficients c(^„(A) are a subset of the short-distance coefficients 
^mn- V'^ s-iid X s-re two-component Pauli spinors that create a heavy quark and a heavy 
antiquark, respectively, is a Pauli matrix, and D is the spatial part of the covariant 
derivative acting to the left and right anti-symmetrically. Note that there is no sum over 
A on the right side of Eq. (15). All of the three- vector quantities in the NRQCD matrix 
elements for the Hi are defined in the Pi rest frame. We will clarify below the meaning of 
the approximation that we have taken to arrive at Eq. (15). 

The short-distance coefficients c^„(A) can be obtained from the full QCD amplitude 
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A'^i^Si, Pi, A; ^So, P2) in Eq. (12). The NRQCD expansion of the full QCD amplitude is 

00 oc 



m=0 n=0 



x(ggC5o)|^n-|^)'"x|o) 

00 00 



= 87Ve£^(gO^(g2) E E Cn(A)qi^"^qr ■ (16) 

m=0 n=0 

In Eq. (16), we use relativistic normalization for the Q and Q states in the computation 
of ^q('^S'i, Pi, A; ^5*0, P2) and in the computation of the NRQCD matrix elements. Conse- 
quently, a factor 4P(gi)P(g2) appears in the second equality of Eq. (16). An additional 
factor 2Nc arises from the spin and color factors of the NRQCD matrix elements. From 
Eq. (16), it is straightforward to calculate the short-distance coefficients c^„(A): 

■Aq^^Si, Pi, A; ^Sq, P2 



Cn(A) = 



1 9™ (9" 



8N,E{q,)E{q2) 



(17) 



Substituting the short-distance coefficients (17) into Eq. (15), one finds that 



00 00 

xEE^-(A)(q^'")j/^(q^"),. 

m=0 n=0 



X 



d 



m=0 n=0 



d 



'^q(^'S'i, Pi, A; ^-S'o, P2) 
AE{q,)E{q2) 



(18) 



qf=q^=0 



Here, the quantities {q )h are ratios of NRQCD matrix elements: 



(J/^(A)|^t(_ll))2n^^.^(A)x|0) 

(J/V'(A)|V'to-.e(A)x|0) 



(r/clV'^xlO) 



(19a) 
(19b) 



We note that (J/V^(A)|V'V . e(A)x|0) and {q 



2m\ 



are independent of the J/%Ij helicity A, 



and there are no sums over A in these quantities. 

Now we can clarify the meaning of the approximation that was taken to arrive at Eq. (15) 
and, consequently, to arrive at Eq. (18). Suppose that we specialize to the Coulomb gauge. 
Then, we can drop the gauge fields in covariant derivatives in the matrix elements in Eq. (18), 
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making errors of relative order v^. The matrix elements are then proportional to derivatives 
of the Coulomb-gauge color-singlet QQ quarkonium wave function at the origin [2]. That 
is, (qf^") is just the 2nth moment of the momentum-space wave function with respect to the 
wave-function momentum (the relative momentum of the Q and Q). Hence, Eq. (18) has the 
interpretation of the convolution of the short-distance amplitude with the momentum-space 
quarkonium wave functions, where the short-distance coefficients have been Taylor expanded 
with respect to the wave- function momenta. Therefore, we see that the approximate NRQCD 
expansion in Eqs. (15) and (18) includes all of the relativistic corrections that are contained 
in the color-singlet QQ quarkonium wave function, up to the ultraviolet cutoff of the NRQCD 
matrix elements.^ 



V. RESUMMATION 

In Ref. [8], a method was presented for resumming a class of relativistic corrections to 
the color-singlet (S'-wave amplitudes that appear in the production and decay of yS-wave 
quarkonium states. The key to the resummation is an expression that relates the (S-wave 
color-singlet matrix elements of higher orders in v to the matrix elements of relative orders 

and v^: 

{q''')H = {qTh- (20) 

The relation (20) is derived in the approximation in which the Q and Q interact only through 
the leading spin-independent potential. Consequently, the relation (20) is accurate up to 
corrections of relative order f ^.^ 

The amplitude (18) is a function of the ratios and (q^'*)??^. Applying the relation 

(20) to Eq. (15), one obtains the resummed expression 

A''[J/i:{X)+Vc] = ^{JmX)\^^cr-e{X)x\0){Vc\^h\0) 



X 



V^V^ 3 1 



• (21) 

qf=(q^)j/v..92=(q^>r,c 



^ We note that, in the case; of dimensionally regulated NRQCD matrix elements, pure power ultraviolet 
divergences in the matrix elements are set to zero. Hence, the effects of pure-power-divergent contributions 
are absent in the resummation. 

^ The derivation involves specializing to the Coulomb gauge and replacing covariant derivatives in operators 
with ordinary derivatives. This approximation also introduces errors of relative order v^. 
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The expression (21) resums those relativistic corrections that are contained in the QQ 
quarkonium wave function in the leading-potential model for the wave function. We note 
that, because the relation (20) is accurate only up to corrections of relative order v^, the use 
of the resummed expression (21) generally does not improve the nominal accuracy over that 
which one would obtain by retaining only corrections through relative order v"^ in Eq. (15). 
The exception to this is the situation in which the short-distance coefficients in Eq. (17) 
grow rapidly with m or n. Then the terms of nominally higher order in v in Eq. (15) can 
have numerical values that are comparable to or larger than the numerical value of the term 
of nominal order v^. In that situation, the resummed expression can give an improved esti- 
mate of the amplitude. The resummed expression may also give an indication of the rate of 
convergence of the v expansion. In any case, it is generally useful to include a well-defined 
set of higher-order contributions in a calculation whenever possible. 



VI. CHOICE OF FRAME AND CO-ORDINATE SYSTEM 

In calculating Aq(^Si, Pi, A; ^Sq, P2), it is convenient to speciahze to the e+e" CM frame, 
to choose a particular coordinate system, and to choose a particular convention for the 
polarization vectors of the QQii^Si) states for the various helicities. We make these choices 
as follows: 

/ci = -^(l, + sin^,0, + cos^), (22a) 

A;2 = ^(l,-sin^,0,-cose), (22b) 

P* = (Pi,0,0,+Pcm), (22c) 

P* = (£;2,0,0,-Pcm), (22d) 

6*(0) = , / (PcM, 0,0,^1), (22e) 

V -^1 ~ ^CM 

e*(±) =T-^(0,l,Ti,0). (22f) 

Here, the angle is the scattering angle, four-vectors are written as v — (v°, v^, v^, v^), and 

_ AV^(5,mf,mi) 
PcM - -z-^ > I2daj 



Ei=\PlM + rhl (23b) 
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where X{x,y,z) — x"^ + y'^ + z'^ — 2{xy + yz + zx). We have used the notation P* to 
distinguish the particular values of these quantities in the e'^e" CM frame from the values 
in the quarkonium rest frame (Sec. Ill A). As with the quantities mi and m2 in Eq. (14), 
rhi and 1712 represent the J/ip and rjc masses, respectively. We will specify below how these 
are chosen for various parts of the calculation. 

Now let us write expressions for the relative momenta qi and q2 in the e+e" CM frame. 
In the quarkonium rest frame, Qi is given by 

Qi = l^il (0, sin 9i cos 4>i, sin 9i sin 0j, cos 9i), (24) 

where 9i and 0^ are the polar and azimuthal angles of qi. Boosting Eq. (24) from the Pi rest 
frame to the e'^e~ CM frame, one obtains 

Qi — \Qi\ i+liPi cos 9i , sin 9i cos 0i , sin 9i sin 0i , 71 cos 9i ) , (25a) 

12 = 1 92 1 (-72/92 cos 612, sin ^2 cos 02, sin 6^2 sin 02, 72 cos 6^2), (25b) 

where 

7i = Ei/^Ef-P^^, (26a) 

lA = Pcu/^Ef-P^^. (26b) 

Note that \qi\ — \f—^ is the magnitude of the three- vector, not in the e+e" CM frame, but 
in the Pi rest frame. 

It follows from Eq. (22) and the analogue of Eq. (3) for ^^(^-Si, Pi, A; ^S^, P2) that 

^^(35i,Pi,0;%,P2) =0, (27a) 
A'i^CSi,Pi,±;'So,P2) = ±AQPcMV^e*'*(±). (27b) 

It is efficient to determine Aq by carrying out the computation of the amplitude 
^^(^5*1, Pi, A; ^5*0, P2) for one value of /i and one value of A such that e*'^(A) is nonzero. 

VII. CROSS SECTION 

Making use of Eq. (3) and the exphcit choices of helicity states in Eq. (22), we find that 

A^[J/tlj{0)+rjc]^0, (28a) 
A''[J/^{±) + Vc] = ±APcMVie*n±)- (28b) 
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Comparing Eq. (28) with Eq. (27) and making use of the resummed NRQCD expansion in 
Eq. (21), we see that 



2E{q,)2E{q2) 



A. 



Q 



(29) 



Cli = {<l^)j/^,<ll={<l^)vc 

A4{X), the S'-matrix element for the process e+e~ — > J/V' + rjc, is defined in Eq. (1). By 
making use of Eqs. (3) and (28), one can evaluate A4{X): 



M{0) = 0, 



(30a) 
(30b) 



The squared helicity amplitudes, summed over the spin states s"*" = ±1/2 and s — ±1/2 
of the e+ and e~, respectively, can be obtained by using Eqs. (22) and (30): 



(31) 



s±=±l/2 



which lead to 



E 1-^(0)1' = 0' 



s±=±l/2 



E \Mi±)\' = ey\A\'PU^ + cos'9). 



(32a) 
(32b) 



s±=±l/2 



Averaging the squared helicity amplitude (32) over the lepton spins, dividing by the flux 
2s, and integrating over the two-body phase space, we obtain the total cross section for 
e+e~ J/ip + Tjc' 

a[e+e-^J/V. + r;.] = lxlx$,|_'^^E E l-^^l' 



A=± s±=±l/2 



An' 



3Nls 



2mi 2m2 



\Ar 



9l = (q^)j/V"92=('J^)'7c 



(33) 



where 



|(J/V(A)|V'V.e(A)x|0)f , 
= \{rjMO)\\ 



(34a) 
(34b) 
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and $2 is the two-body phase space 

^'"8^^'^'^^''"'/^''"^=^- ^^^^ 
Note that we use the physical masses for the J/ if) and the r]c in the phase space (35). 

VIII. VMD TREATMENT OF THE PHOTON-FRAGMENTATION AMPLITUDE 



In the part of the amphtude that comes from the photon-fragmentation diagrams of the 
type in Fig. 2, we can reduce the theoretical uncertainty by making use of the VMD method 
to calculate the fragmentation of the 7* into J/ip (Ref. [13]). In Ref. [13], the process 
Y ^ J/il) has been calculated using the VMD method. Using Eq. (3) of Ref. [13], we find 
that the 7* — > J/'^ couphng is 

^./.= (^lJ^r[W-^^r]j . (36) 

In order to implement the VMD calculation we must make the following substitutions in 
the NRQCD calculation of the photon-fragmentation diagrams: 

ec^j2mi{Oi)j/^ gj/^, (37a) 
'IV{(7^® l)[n3(]9i,pi,A)®7ri]} 



^/2Nl2E{q^) 

where gj/^ is defined in Eq. (36). 



<(A), (37b) 



IX. CHOICE OF THE QUARKONIUM MASSES 

We now specify our choices of the quarkonium masses in our computation. In computing 
diagrams involving on-shell quarks, such as those in Fig. 1, it is generally necessary, in order 
to maintain gauge invariance, to choose the quarkonium masses so as to respect the on-shell 
condition. Hence, we generally must choose rfij = 2E(qi) = 2 ml + qf in Eq. (23) in 
working out the kinematics. 

An exception to this is in the computation of the photon-fragmentation diagrams of the 
type in Fig. 2. In this case, if we make use of the VMD method for calculating the amphtude. 
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we can maintain gauge invariance even if we take rhi to be rrij/^, the physical J/ip mass. 
We still must choose m2 = 2£'(g'2) for the rjc mass, however. Since one generally reduces 
theoretical uncertainties by eliminating 2mc in favor of mj/^, we choose mi = rrij/^ in the 
VMD calculation of the photon-fragmentation diagrams. 

The factor \l2m\\l2m2 in Eq. (21) arises from the relativistic normalizations of the states. 
In this case, we choose = 2E{qi). It turns out that this choice leads to a near cancellation 
of the dependence on rric in the amplitude at leading order in v (Ref. [3]). Thus, this choice 
reduces the theoretical uncertainties that arise from the uncertainty in rric- 

As we have already noted, we use the physical quarkonium masses in computing the phase 
space in Eq. (35). At first sight, this choice might appear to be inconsistent with the choice 
rrii = 2E[qi) in Eq. (21), since the factors rrii in Eq. (21) arise from the normalizations of 
the states, which also enter into the phase space. The choices that we have made amount 
to multiplying the amplitude by the factors ■^y2E(qi)/mj/^ and ^/2E{^^)Jm^^. At the level 
of precision in v to which we work, these factors are equivalent to unity. 

X. INTERFERENCE WITH THE NLO AMPLITUDE 

As we have mentioned, the corrections to a[e^e^ J /ip + ?7c] at NLO in have been 
calculated in Ref. [7]. Because the amplitude for the relativistic corrections has the same 
phase as the amplitude at leading order in i>, we can infer, from the results of Ref. [7], the 
contribution to the cross section of the interference between the amplitude at NLO in ag 
and the amplitude for the relativistic corrections. 

First, let us define some notation. When we discuss cross sections a and reduced hadronic 
amplitudes A [Eq. (3)], a subscript indicates that the quantity is computed at leading order 
in f , a subscript v indicates that the quantity is resummed to all orders in a subscript 
NLO on A indicates the contribution to A at NLO in as, and a subscript NLO on a indicates 
the sum of the contributions to a through NLO in CKg. A superscript QCD indicates that 
only QCD contributions to the hadronic amplitude have been included. The absence of a 
superscript QCD indicates that both QCD and pure QED contributions to the hadronic 
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amplitude have been included. Using this notation, we have 



(70 = N\Ao\\ 



QCD 



(38a) 
(38b) 
(38c) 
(38d) 



where the normalization factor M is that of Eq. (33) and is defined by 



3s 



-ec«'^CM*2, 



(39) 



with PcM defined in Eq. (23) and $2 defined in Eq. (35). 

The quantity ct^nlo computed in Ref . [7] . On the other hand, we wish to compute the 
quantity 

(7tot = ^ 2Re(A^;;S)) ■ 

Using the fact that and jQ^^ have the same phase, we can write 



(40) 



ORpM 4*QCD\ _ pp / .QCD .*QCD \ 

^0 

QCD _ QCD 
_ -I- y^"0,NLO "0 

M ^ I QCD 



(41) 



Thus, (Ttot can be expressed in terms of cr^, (7^*"°, and c^nlo- 



QCD 



QCD 



(T'tot = + 



^0, NLO ^0 



QCD 



QCD 



(42) 



XI. RESULTS 



In this section, we present our numerical results. 

We compute Aq in Eq. (27) from the Feynman diagrams in Figs. 1 and 2, making use of 
the spin and color projectors, as described in Section III. We then carry out the projection 
onto the (S-wave states by performing the integration over the angles of q\ and q-i numerically, 
as indicated in Eq. (13). We substitute Aq into Eq. (29) to obtain A and substitute A into 
Eq. (33) to obtain the cross section cr„, which includes the resummed relativistic corrections. 
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We compute ctq, the cross section at leading order in v, by setting — — in the 
expressions for ay. 

In carrying out this calculation, we make use of the matrix elements {Oi)j/^ and (Ci),,^ 
and the ratios of matrix elements {q'^)j/ip and {q^)r,c from Tables 1 and 111 of Ref. [11]. In 
Ref. [11], various uncertainties were associated with these quantities. The uncertainties are 
correlated to varying degrees among the quantities. We recount the uncertainties here. 

There are theoretical uncertainties in the values of {q^)j/.^ and {q^)nc that arise from 
the fact that the leading-potential approximation that is used in Ref. [11] is accurate only 
up to corrections of relative order v^. These uncertainties are denoted by A{q'^)j/,^, and 
A(qr^)^^, respectively. There arc uncertainties that arise from the scale uncertainties in 
as and from neglecting next-to-next-to- leading- order (NNLO) corrections to the J/ip and 
rjc electromagnetic widths. They are denoted by ANNLOj/^ and ANNLO^^, respectively. 
There are also uncertainties that are associated with the heavy-quark mass rric, the string 
tension a, and the uncertainties in the experimental measurements of r[J/ip — > e'^e"] and 
r[?7c — 77]- These uncertainties are denoted by Arric, Act, AFj/^, and AF^^, respectively. 
Finally, there is an uncertainty that is associated with the use of the heavy-quark spin 
symmetry to combine the values of the rjc matrix elements that were obtained from T[r)c 
77] with those that were obtained from T[J/ip e+e"]. It is denoted by Av"^. 

The uncertainty estimates in Ref. [11] make use of the standard NRQCD power-counting 
(velocity-scaling) rules [2]. Alternative power-counting rules have been proposed [14-16], and 
the use of these rules would lead to estimates for (q^) j/^i and (q^),,^ that are of relative order 
unity. However, lattice calculations [17-20] support the notion that the standard NRQCD 
power-counting rules give an upper bound on the uncertainties. Therefore, we make use of 
the uncertainty estimates of Ref. [11]. (See Ref. [11] for a more detailed discussion of these 
issues.) 

The calculation requires some additional inputs. We take y/s — 10.58 GeV. In order to 
maintain consistency with the calculation at NLO in as in Ref. [7], we take rric to be the 
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one-loop pole mass. The specific numerical value that we use is^ 

rric = 1.4 ± 0.2 GeV. (43) 

This choice of numerical value corresponds to the one in Ref. [11], and so, in determining 
the uncertainties that arise from the uncertainty in rric, we are able to make use of the 
dependences of the matrix elements and ratios of matrix elements on rric that are computed 
in Ref. [11]. For the electronic width of the J/ip, which enters into the calculation of the 
VMD couphng gj/^ from Eq. (36), we take 

rlJ/ijj e+e-] = 5.55 ± 0.14 ± 0.02 keV. (44) 

For the strong and electromagnetic couphngs we take^ 



and 



q;,(10.58/4 GeV) = 0.26, 
q;,(10.58/2 GeV) = 0.21, 
q;,(10.58 GeV) = 0.17, 

q;(10.58/4 GeV) = (132.9)-\ 
q;(10.58/2 GeV) = (131.9)-\ 
a(10.58 GeV) = (130.9)-\ 
a{mj/^) = (132.6)-\ 



(45a) 
(45b) 
(45c) 

(46a) 
(46b) 
(46c) 
(46d) 



We determine the central value of the scale for each coupling from the momentum transfer 
at the relevant vertex. Let us call the virtual photon that connects the lepton and c-quark 
lines photon 1, the non-fragmentation virtual photon that connects c-quark lines in Fig. 1 
photon 2, and the virtual photon in the fragmentation diagrams in Fig. 2 photon 3. At 
virtual-gluon vertices and at photon-2 vertices, we take the scale to be half the CM energy; 

* The most recent compilation of the Particle Data Group [21] suggests that the uncertainty in rUc may be 
a factor of two smaller than the uncertainty that we have used here. However, since it is not clear that 
the systematic errors are well understood in the various determinations that enter into that compilation, 
we make a conservative choice of error bars here. 

^ We compute and a at each scale by making use of the code GLOBAL ANALYSIS OF PARTICLE 
PROPERTIES (GAPP) [22]. 
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TABLE I: The cross sections in units of fb that were obtained by Zhang, Gao, and Chao [7, 23]. The 
first row below the headings contains the cross sections for central values of rric and /i. Subsequent 
rows contain the cross sections for the plus and minus variations of rric and fx with respect to 
their uncertainties. The strong coupling constant was taken to be a^*^*-" (10.6/4 GeV) = 0.273, 
a^GC(io.6/2 GeV) = 0.211, and a^^^{10.6 GeV) = 0.174 (Ref. [23]). 

PflSP n-ZGC -.ZGC 

i^ase (7q ctq^ ]^lo 

central 5.8 12.6 

+Amc 4.8 9.7 

-Arric 6.7 15.7 

+A// 3.9 8.9 

-An 9.6 19.6 

at photon- 1 vertices, we take the scale to be the CM energy; at photon-3 vertices, we take 
the scale to be the J/ip mass. In calculating the VMD couphng gj/^ from Eq. (36), we also 
take the scale of the virtual- photon vertices to be the J/i/j mass. These choices are consistent 
with those in Ref. [11]. We also use mj/^ = 3.096916 GeV and m^^ = 2.9798 GeV [21]. 

In computing the cross section (Jtot, which includes relativistic corrections, corrections of 
NLO in as, and the interference between them, we make use of Eq. (42). We compute the 
quantity (ct^nlg ~ ^o'^^)/ \J ^^^^ by making use of results [23] that have been provided by 
the authors of Ref. [7]. These results are shown in Table I. Here we have introduced an ad- 
ditional uncertainty A/i, which accounts for effects of the uncertainty in the renormalization 
scale in the NLO calculation. This uncertainty is determined by varying the renormalization 
scale by a factor of two above and below its central value of 10.58/2 GcV. In making use of 
the results from the authors of Ref. [7], we rescale the values of q;^, a, the matrix elements 
{Oijj/^ and phase space so that they conform to our choices. The rescahng 

is carried out as follows: 

QCD _^QCD ZGC ^ZGC 

-^Tp T^. , (47) 



where the superscript ZGC indicates the value that was given by the authors of Ref. [7] 
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The scaling factor p is defined by 

a(10.58 GeV) = 1/130.9\ Y V $2 



(48) 



The values for and (Oi) are given by 

^f^°-^x^/l-(^) , (49a) 
3 




{Oi)'ff? = 7^ X 0.978 GeV^ = 0.467 GeV^ (49b) 



Our numerical results are shown in Table II. The first row below the headings gives 
the central values of the matrix elements, ratios of matrix elements, and the cross sections. 
Subsequent rows contain the maximum and minimum values of each of these quantities that 
are obtained by varying the input parameters with respect to each of the uncertainties that 
we have described. The matrix elements and the ratios, as well as their variations with 
respect to each uncertainty, arc taken from Tables I and III of Rcf. [11]. The deviations 
from the central values, given in the same order as the rows in Table II, are as follows: 

rr - (\ 4+0.1+0.5+0.1+1.5+0.3+1.1+0.4+0.5+0.7 fU _ c 4+2.I r- (^-n. \ 
(70 — D.4_o^_o 5_o^_^^_o 2_i,i_o.3-0.4-0.8 " 9 ID, 

rr — q+0.5+2.5+0.4+2.0+0.3+1.5+0.9+0.7+1.0 fu _ n q+3.9 r< ( Kf\U\ 

Uy — y-<J_o.5-1.7-0.4-1.5-0.3-1.5-0.9-0.6-l.l — ^■•^-3.2 '-'^1 yo\JD) 

rr —17 c+0.8+5.3+0.7+3.9+0.7+2.8+1.6+1.4+1.9 ru, _ 1 7 fi+7.8 ru, (rr^^\ 
C'tot — J- '■O_Q 9_3 7_0 7_3 0-0.7-2.9-1.5-1.1-2.0 ~ -'-'■"-6.3 \0\J(.) 

In the result for crtot above, we have not included the uncertainty A// that arises from varying 
the renormahzation scale. That uncertainty is igig fb. This is, perhaps, an overestimate of 
the uncertainty from uncalculated corrections of higher order in as and v^, since it assumes 
that our choice of renormahzation scale may be wrong by as much as a factor of two. Alter- 
natively, one could estimate the uncertainty that arises from uncalculated corrections in the 
following way. One could take for the uncertainty associated with uncalculated corrections 
of NNLO in a, to be the quantity ANNLO = Q;,((7Jnlo - <^o^^°) ~ 1-32 fb, and one could 
take for the uncertainty associated with uncalculated corrections of NLO in ag and NLO in 
v"^ the quantity ANLO-'U^ = '^^(f''o!NLO ~ ^?*^°) ~ fb- (Uncertainties of relative order v'^ 
are already included in A{q^)j/^ and A(qr^)^^.) If we add the uncertainty A/x in quadrature 
with the other uncertainties, then we obtain 

atot = n.Qll'i" fb. (51) 
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TABLE II: The matrix elements (C'i)j/^ and {Oi)n^ in units of GeV^, the ratios of matrix elements 
{q^)j/^ and {q^)nc in units of GeV^, and the cross sections ctq, cr^, and atot in units of fb. The first 
row below the headings contains central values for the matrix elements, the ratios, and the cross 
sections. Subsequent rows contain the maximum and minimum values of each of these quantities 
that are obtained by varying the input parameters with respect to each uncertainty. 



Case 












(7„ 


O"tot 


central 


0.440 


0.441 


0.437 


0.442 


6.4 


9.3 


17.6 




0.450 


0.573 


0.437 


0.442 


6.5 


9.8 


18.4 




0.430 


0.308 


0.437 


0.442 


6.3 


8.8 


16.7 


+Amc 


0.433 


0.443 


0.470 


0.430 


6.0 


7.6 


13.9 


— Amc 


0.451 


0.437 


0.413 


0.450 


6.9 


11.8 


22.8 


+A(7 


0.443 


0.482 


0.444 


0.482 


6.6 


9.7 


18.3 


-A(7 


0.437 


0.400 


0.431 


0.403 


6.3 


8.9 


16.9 


+ANNLOj/^ 


0.504 


0.419 


0.473 


0.429 


7.9 


11.3 


21.5 


-ANNLOj/^ 


0.387 


0.459 


0.408 


0.452 


5.3 


7.8 


14.6 




0.451 


0.437 


0.443 


0.440 


6.7 


9.6 


18.2 


-AFj/^ 


0.429 


0.444 


0.431 


0.444 


6.2 


9.0 


16.9 


+Au2 


0.440 


0.441 


0.511 


0.417 


7.5 


10.8 


20.4 


-Ai;2 


0.440 


0.441 


0.364 


0.467 


5.3 


7.8 


14.7 


+A(g2),. 


0.440 


0.441 


0.461 


0.574 


6.8 


10.2 


19.1 


-new 


0.440 


0.441 


0.414 


0.309 


6.1 


8.4 


16.1 


+ANNLO^, 


0.440 


0.441 


0.474 


0.429 


7.0 


10.0 


19.0 


-ANNLO^^ 


0.440 


0.441 


0.408 


0.452 


6.0 


8.7 


16.4 


+Ar^e 


0.440 


0.441 


0.487 


0.425 


7.2 


10.3 


19.5 


-AT,, 


0.440 


0.441 


0.385 


0.460 


5.6 


8.2 


15.5 




0.440 


0.441 


0.437 


0.442 


4.4 


6.3 


12.3 


-All 


0.440 


0.441 


0.437 


0.442 


9.5 


13.9 


25.0 
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On the other hand, if we add ANNLO and ANLO-^;^ in quadrature with the other uncer- 
tainties, then we obtain 

(Ttot = 17.611^ fb. (52) 

In addition to the uncertainties that we have included in Eqs. (51) and (52), there are 
uncertainties that are associated with the NRQCD factorization formula. A rigorous proof 
of NRQCD factorization for crfe+e" — > J/t/j + rjc] does not exist. However, it seems likely, on 
the basis of existing work on proving NRQCD factorization for other production processes 
[24], that the corrections to the factorization formula are of order m^/(s/4) ^ 34%, where 
TTiH is the mass of either of the heavy quarkonia. 

The various contributions to atot are as follows. The cross section at leading order in 
as and v, ctq, contributes about 6.4 fb, of which about 1.0 fb comes from the pure QED 
corrections. The direct relativistic corrections that are associated with the process e'^e' — > 
J/t/j + ric contribute about 2.9 fb. The corrections of NLO in contribute about 6.9 fb, 
including the interference with the pure QED contribution. The interference between the 
relativistic corrections and the corrections of NLO in contributes about 1.4 fb. 

We have examined our numerical calculation in the limits {q'^)j/^ and {q^).,,^ — > and 
find that it agrees with the analytic results in Refs. [3] and [10] for the order-v^ corrections 
to a[e'^e'' — > J/ip + r]c]- 

The direct relativistic corrections to the process e'^e~ — > J/t/j+rjc itself are modest in size. 
ay is about 45% larger than uq, but uo already contains an implicit relativistic correction 
factor of 0.96 that arises from the use of the hadron masses, rather than 2mc, in the phase 
space. Hence, the enhancement from the direct relativistic correction is about 40%. If we 
use the hadron masses in the phase space and keep only the order- relativistic corrections 
to the squared amplitude, then we find that the direct relativistic corrections increase the 
cross section by about 45%. Thus, we see that the effects of resummation are not large, 
suggesting that the velocity expansion of NRQCD is converging well in this case. 

As wc have mentioned previously, the rcsummcd result contains all of the corrections 
that are associated with the momentum-space QQ quarkonium wave function in the leading- 
potential approximation, up to the ultraviolet cutoff of the NRQCD matrix elements. Hence, 
the modest size of the relativistic corrections supports the conclusion in Ref. [25] that the 
effects of the finite width of the momentum-space QQ wave function are not dramatic, once 
one excludes contributions from the large-momentum tails of the wave function. Those 
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contributions are included in the NRQCD formalism in the corrections of higher order in 

OLs- 

The use of the VMD method, rather than the NRQCD method, in calculating the frag- 
mentation amplitude in the pure QED contribution, has a small effect on the central value 
of the cross section. The use of the VMD method shifts the central value of cr^ down by 
about 3%. 

XII. COMPARISON WITH PREVIOUS CALCULATIONS 

Let us now compare our results with some of those from previous calculations. 

As we have already mentioned, the contribution to (j[e"^e~ J/ip + rj^ at leading order in 
as and v was first calculated in Refs. [3] and [4]. There are some differences in these results, 
owing to different choices of input parameters and the inclusion of pure QED corrections 
in Ref. [3]. Let us focus on Ref. [3], since the calculation in that paper is closer to the 
present one in terms of input parameters and the treatment of pure QED corrections. The 
result in Ref. [3] is a[e^e~ — > J/t/j + rjc] — 3.78 ± 1.26 fb. This result should be compared 
with our result for uo, which is about 70% larger. This difference arises essentially because 
we have used the values for {Oi)j/^ and {Oi)r,^ from Ref. [11] (see Table 11), while the 
authors of Ref. [3] have used (Ci) j/^ = {Oi)^^^ = 0.335 GeV^. This substantial difference in 
the values of the matrix elements arises largely from the inclusion of relativistic corrections 
to the electromagnetic decay widths of the J/t/j and the rjc in analyses of Ref. [11]. The 
relativistic corrections increase the sizes of and (Ci)r;c by about 31%, where we are 

comparing in both instances with the matrix element that is extracted from T[J/iIj e~^e~] 
in Ref. [3]. The changes in the values of the matrix elements lead to a 72% change in the 
cross section. Other small differences in our calculation relative to that in Ref. [3] arise 
from using the VMD method to calculate the fragmentation contribution to the pure QED 
amplitude (about —8%), from the use of the physical masses for the J/ip and the rjc in 
the phase space (about —4%), and from taking into account the effects of the running of a 
(about 10%). The error bar in the result of Ref. [3] takes into account only the uncertainty 
Arric- It is more than twice the size of the Arric error bar in (Tq in our calculation. The error 
bar in our calculation is reduced because the Arric uncertainty in the matrix elements in 
Ref. [11] was reduced by replacing certain factors of 2mc with mj/^. 
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In Ref. [9], a result crtot = 17.5 ± 5.7 fb is given. Our calculation contains a number of 
refinements in comparison with that of Ref. [9]. Among them are the use of the improved 
results for the matrix elements in Ref. [11], the use of the VMD method to calculate the 
fragmentation contribution to the pure QED amplitude, the inclusion of the effects of the 
running of a, and the precise calculation of the interference between the relativistic correc- 
tions and the corrections of NLO in cts, rather than the use of an overall K factor to account 
for the corrections of NLO in ctg. The effects of these refinements cancel almost exactly in 
the central value for the cross section. The error bars in the result of Ref. [9] include only 
the uncertainties Amc, A(g^)j/^, and /S.{q^)ri^ and are, therefore, somewhat smaller than 
the error bars that wc report here. 

Wc can also compare our results with those of Ref. [10]. In that work, the quantities 
{Oi)nc^ and {q'^)j/,p{Oi)j/^ = (q^),,,(Ci)r?c were determined by comparing theoret- 
ical expressions for r[J/t/j — > e+e~], r[r]c — > 77], and T[J/t/; — > hght hadrons] with the ex- 
perimental measurements of those widths. The resulting values are {Oi)^J^ — 0.573 GeV^, 
{Oi)^J^ = 0.432 GeV^ {q^)^J^ = 0.202 GeV^ and {q^)^J^ = 0.268 GeV^ where the su- 
perscript HFC denotes the value that was given in Ref. [10]. The value of {Oi)j/^ is about 
30% larger than the one that we employ, and the value of {Oi)^^ is about 1% smaller. The 
values of and {q'^)r]c are smaller than the values that we use by about 54% and 39%, 

respectively, and are considerably smaller than expectations from the NRQCD velocity- 
scahng rules. As was discussed in Ref. [11], the smaller values of and (q^)r,^ arise 

in Ref. [10] because the theoretical expression for V[J/ip — > hght hadrons] contains a very 
large relativistic correction. We regard this as an indication that the velocity-expansion for 
that process is not under control. 

The authors of Ref. [10] find that the direct relativistic corrections to (T[e~^e~ — * J/ip + rjc] 
enhance the cross section by about 26%. If we ignore the effects of pure QED contributions 
and resummation, we find an enhancement from direct relativistic corrections of about 56%. 
The difference presumably arises from the use of smaller values of and (q^),,^ in 

Ref. [10]. 

In Ref. [10], the central value for the total cross section is o-^f^ = 20.04 fb. This result 
does not include the pure QED contribution, the contribution from the interference between 
the corrections of NLO in ag and the relativistic corrections, and the effects of resummation. 
The corresponding quantity in our calculation is 14.7 fb. Thus, we see that the result of 
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Ref. [10] is 37% larger than ours. The main sources of this difference are the use of a larger 
value of which would increase our cross section by about 30%, the use of a larger 

value of the strong coupling (a^^^ = 0.2592), which would increase our cross by about 
47%,^ the use of a smaller value of the electromagnetic coupling (a^^^ = 1/137), which 
would decrease our cross section by about 9%, the use of a larger value of the charm-quark 
mass (m^^*" — 1.5 GeV), which would decrease our cross section at fixed values of the 
NRQCD matrix elements by about 12%, and the use of smaller values of (q^) j/^ and {q'^)r)c: 
which would decrease our cross section by about 9%. 

In Ref. [10], the dependence of the cross section on rric is given. As rric is varied from 
1.4 GcV to 1.6 GeV, a change in the cross section of +37% is found. In contrast, for this 
variation in tjIc, we find a change in the cross section of —30%. Presumably the difference 
arises because, in the method that is used to determine {Oi)j/^ and (Ci)?,^ in Ref. [10], 
those quantities are proportional to m^. In the method that is used in Ref. [11] to determine 
{Oi)j/^ and {Oi)r,^, the dependence of those quantities on rric is much milder, partly because 
some factors of 2mc are replaced with rrij/^ in the theoretical expressions. The authors of 
Ref. [10] have not estimated the sizes of uncertainties that arise from other sources, and so 
it is not clear whether their method of calculation leads to a more precise prediction for the 
cross section than the one that we have used. 



XIII. SUMMARY AND DISCUSSION 



For a number of years, the discrepancy between theoretical predictions for the exclusive 
double-charmonium cross section a[e'^e~ J/ip + rjc] and experimental measurements has 
posed a significant challenge to our understanding of quarkonium production. Changes in 
the measured values of the cross section have reduced the discrepancy somewhat [5, 6]. 
More recently, calculations of the corrections of NLO in ag [7] and relativistic corrections 
[9, 10] have increased the theoretical prediction for the cross section by almost an order 
of magnitude. The shifts in the theoretical and experimental central values for the cross 

^ In Rof. [10] the cross section at NLO in was computed using a,, = 0.2592 and fj, = 3.00 GeV (Ref. [23]). 
In order to find the effect of this choice of ag and /i on our calculation, we compute (Tq'~''~^ and ctq^^J^q at 
Qg = 0.2592 and fi = 3.00, using the cross sections in the last row of Table I as inputs. We then use these 
values for cto'^^ and ctq^^q to evaluate Eq. (47). 
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section have resolved the outstanding discrepancy. However, in the absence of an analysis 
of the theoretical uncertainties, the meaning of the apparent agreement between theory and 
experiment is unclear. 

In this paper, we have carried out a new computation of the relativistic corrections to 
a[e~^e~ — >■ J/ip + rjd, with the goals of adding certain refinements to the calculation and 
making a more precise estimate of the theoretical uncertainties. Some of the refinements, 
relative to the calculation of Ref. [9], are the use of the VMD method to calculate the 
fragmentation contribution to the pure QED amplitude, the inclusion of the effects of the 
running of a, and the inclusion of a precise calculation of the interference between the 
relativistic corrections and the corrections of NLO in a^, as opposed to the use of a simple K- 
factor estimate. A further significant refinement in our calculation is the use of an improved 
determination of the relevant NRQCD matrix elements at leading order in v"^ and at NLO in 
v'^ (Ref. [11]). This determination includes an analysis of the correlated uncertainties in the 
matrix elements. Our calculation exploits this information to give a much more complete 
estimate of the uncertainties than was given in Ref. [9] . 

Our calculation differs from the one in Ref. [10] in that we include pure QED corrections, 
we take into account the effects of the running of a, we include the interference between the 
relativistic corrections and the corrections of NLO in ag, and we resum a class of relativistic 
corrections. As we discuss in Section XII, the calculation of Ref. [10] makes use of matrix 
elements that differ significantly in numerical value from those that we use. Ref. [10] includes 
a discussion of scale uncertainties and the effect of the uncertainty in rric, but does not provide 
an overall error bar for the cross section. 

In our calculation, the relativistic corrections to a[e'^e~ ^ J/ip + r]c\ arise from two 
sources. The first, direct source consists of the relativistic corrections to the process e+e" — > 
J/i/j + r]c itself. These increase the cross section by about 40%. The second, indirect source 
of relativistic corrections derives from the relativistic corrections to the electromagnetic 
decay widths of the J/ip and the rjc, which enter into the matrix-element determinations of 
Ref. [11]. These corrections increase the cross section by about 88% (Ref. [11]). (Other, 
smaller corrections result in a net change in (Oi) j/^{Oi)rj^ of 72% relative to the value of 
(Oi) j/ip{Oi)rj^, that was used in the calculation of Ref. [3].) The inclusion of corrections of 
NLO in as further increases the cross section by about 89%, of which about 15% comes from 
the interference between the relativistic corrections and the corrections of NLO in q;^. 
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Our principal results are given in Eqs. (51) and (52). In the former result, the uncer- 
tainties that arise from uncalculated higher-order corrections are estimated by varying the 
renormalization scale. In the latter result, those uncertainties are assumed to be given by 
their nominal sizes, namely, Us and f ^ times the contribution to the cross section of NLO in 
as- In addition, there are uncertainties that result from the use of the NRQCD factorization 
formula for the cross section, which we estimate to be about 34%. 

The central value for a[e^e~ — > J/V' + 77c] that we obtain is essentially the same as that 
of Ref. [9]. The effects of the various refinements that we have mentioned largely cancel. 
However, some of the refinements allow us to constrain the theoretical uncertainties more 
tightly. Because we have included more sources of uncertainty in our estimates, our error 
bars are significantly larger than those in Ref. [9]. 

Our results for crfe+e" J/ip + rjc] also agree, within uncertainties, with the result of 
Ref. [10]. To some extent, the effects of our use of different values of the matrix elements and 
other input parameters are canceled by our inclusion of additional corrections. In Ref. [10] , 
the dependence of the cross section on rric is given. That dependence is similar in magnitude 
but opposite in sign to the one that we find, presumably because the authors of Ref. [10] use 
a method to determine the NRQCD matrix elements that is quite different from the method 
in Ref. [11]. The authors of Ref. [10] have not estimated other uncertainties, and so it is not 
clear whether their method of calculation yields a result that is more precise or less precise 
than ours. 

As we have mentioned, in our calculation, we resum a class of relativistic corrections 
to all orders in v. These corrections include all of the relativistic corrections that are 
contained in the color-singlet QQ quarkonium wave function, up to the ultraviolet cutoff of 
the NRQCD matrix elements. The effect of the resummation beyond relative order f ^ is 
small, indicating that the velocity expansion converges well for this process. The fact that 
the direct relativistic corrections are modest in size supports the conclusion in Ref. [25] that 
the effects of the finite width of the momentum-space QQ wave function are not dramatic, 
once one excludes contributions from the large-momentum tails of the wave function that 
are contained in corrections of higher order in a^- 

Let us discuss the prospects for decreasing the uncertainties in our calculation. The 
largest uncertainty arises from the uncalculated terms of relative order agv'^ and relative 
order a^. This uncertainty may be as large as t.'^^- A complete calculation of the order- 
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Ugv'^ corrections, which are the larger ones, seems quite feasible. The calculation of the 
corrections of order-a^ would be a major undertaking, but is not out of the question. The 
next largest source of uncertainty arises from the use of the NRQCD factorization formalism 
itself, which may lead to an uncertainty of about 34%. A more thorough understanding of 
the issues that are involved in constructing a rigorous proof of a factorization theorem 
for crfe+e" J/ip + r]c] may lead to a different estimate of these uncertainties. It is also 
conceivable that one could prove a "higher-twist" factorization theorem that would allow one 
to carry out a systematic computation of corrections to the existing NRQCD factorization 
formula. The uncertainties that arise from the use of the NRQCD factorization formalism 
presumably would decrease as the CM energy of the process e+e~ J/ip + rjc increases. 
However, there arc no prospects for measuring a[e'^e~ — >■ J/ip + rjd at higher energies 
in the immediate future. The uncertainty in rric is the next most important source of 
theoretical uncertainty. We estimate the resulting uncertainty in the cross section to be 
-21%- '^^^ expect to see some progress in reducing this uncertainty, particularly from 
lattice determinations of rric- 

Our result for a[e~^e~ — > J/if) + r^c] agrees, within errors, with the measurements of 
the Belle and BABAR experiments. The uncertainties in our result are quite large, and, 
of course, it would be desirable to reduce these uncertainties, so as to sharpen this test 
of the NRQCD factorization approach to quarkonium production. Nevertheless, it seems 
fair to conclude that the long-standing discrepancy between the theoretical prediction for 
(7[e^e~ J/%1) -\- 77c] and the experimental measurements has been resolved. 
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